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Abstract. 

Measurements of correlation functions and their comparison with theoretical models are often employed in 
natural sciences, including astrophysics and cosmology, to determine best-fitting model parameters and their 
confidence regions. Due to a lack of better descriptions, the likelihood function of the correlation function is often 
assumed to be a multi-variate Gaussian. 

Using different methods, we show that correlation functions have to satisfy contraint relations, owing to 
the non-negativity of the power spectrum of the underlying random process. Specifically, for any statistically 
homogeneous and (for more than one spatial dimension) isotropic random field with correlation function (,{x), 
we derive inequalities for the correlation coefficients r„ = ^{nx) / (,{0) (for integer n) of the form r„i < r„ < r„u, 
where the lower and upper bounds on r„ depend on the rj, with j < n, or more explicitly 

H„_ {e(0),5(x),C(2a;), . . . ,e([n - l]x)} < ^{nx) < E„+ {^(0), e(a;), e(2a;), . . .,?([" - l]x)} . 

Explicit expressions for the bounds are obtained for arbitrary n. We show that these constraint equations very 
significantly limit the set of possible correlation functions. For one particular example of a fiducial cosmic shear 
survey, we show that the Gaussian likelihood ellipsoid has a significant spill-over into the region of correlation 
functions forbidden by the aforementioned constraints, rendering the resulting best-fitting model parameters and 
their error region questionable, and indicating the need for a better description of the likelihood function. 

We conduct some simple numerical experiments which explicitly demonstrate the failure of a Gaussian de- 
scription for the likelihood of ^. Instead, the shape of the likelihood function of the correlation coefficients appears 
to follow approximately that of the shape of the bounds on the r„, even if the Gaussian ellipsoid lies well within 
the allowed region. Therefore, we define a non- linear and coupled transformation of the r„, based on these bounds. 
Some numerical experiments then indicate that a Gaussian is a much better description of the likelihood in these 
transformed variables than of the original correlation coefficients - in particular, the full probability distribution 
then lies explicitly in the allowed region. 

For more than one spatial dimension of the random field, the explicit expressions of the bounds on the r„ are 
not optimal. We outline a geometrical method how tighter bounds may be obtained in principle. We illustrate 
this method for a few simple cases; a more general treatment awaits future work. 

Key words, cosmology - gravitational lensing - large-scale structure of the Universe - galaxies: evolution - galaxies: 
statistics 

1. Introduction 

One of the standard ways to obtain constraints on model parameters of a stochastic process is the determination of its 
two-point correlation function ^(a;) from observational data, where x is the separation vector between pairs of points. 
This observed correlation function is then compared with the corresponding correlation function ^{x]p) from a model, 
where p denotes the model parameter (s). A commonly used method for this comparison is the consideration of the 
likelihood function C{^\p), which yields the probability for observing the correlation function ^{x) for a given set of 
parameters p. It is common (see Seljak & Bertschinger 1993 for an application to microwave background anisotropies, 
Fu et al. 2008 for a cosmic shear analysis, or Okumura et al. 2008 for an application to the spatial correlation function 
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of galaxies) to approximate this likelihood by a Gaussian, 



m^ixi)}\p) oc exp 



1 ^ 



(1) 



where it has been assumed that the random field is homogeneous and isotropic, so that the correlation function depends 
only on the absolute value of the separation vector. Furthermore, it has been assumed that the correlation function is 
obtained at discrete points xf, for an actual measurement, one usually has to bin the separation of pairs of points, in 
which case Xi is the central value of the bin. In ([T]), Cov is the covariancc matrix of the correlation function between 
any pair of separations Xi , Xj . 

In a recent paper (Hartlap et al. 2009) we have investigated the likelihood function for the cosmic shear correlation 
function and found that it deviates significantly from a Gaussian. This study relied on numerical ray-tracing simulations 
through the density field obtained from N-body simulations of the large-scale structure in the Universe. 

In this paper, we will show that the likelihood function of the correlation function cannot be a Gaussian. In 
particular, we show that any correlation function obeys strict constraints, which can be expressed as 

S„_ {^(0), C(a^), C(2a:), . . . , ^([n ~ l]x)} < ^(nx) < S„+ {^(0), ^(a;), ^(2a;), . . . , - l]x)} (2) 

for arbitrary x and integer n; these constraints can be derived by several different methods. With one of these methods, 
one can derive explicit equations for the upper and lower bounds in ([2]) for arbitrary values of n. The basic reason for 
the occurrence of such constraints is the non-negativity of the power spectrum, or equivalently, the fact that covariance 
matrices of the values of a random fields at different positions are positive (semi-)definite. 

The outline of the paper is as follows: In Sect.[2l we obtain bounds on the correlation function using the Cauchy- 
Schwarz inequality, as well as making use of the positive definiteness of the covariance matrix of random fields. It turns 
out that the latter method gives tighter constraints on ^; in fact, these constraints are optimal for one-dimensional 
random fields. We show in Sect. [3] that these bounds significantly constrain the set of functions which can possibly 
be correlation function. Whereas the bounds obtained in Sect. [5] are valid for any dimension of the random field, 
they are not optimal in more than one dimension; we consider generalizations to higher-order random fields, and to 
arbitrary combinations of separations Xi in Sect. [4] In Sect. [5] we introduce a non-linear coupled transformation of the 
correlation coefficients based on the bounds; for the case of a one-dimensional field, all combinations of values in these 
transformed quantities correspond to allowed correlation functions (meaning that one can find a non-negative power 
spectrum yielding the corresponding correlations). Hence, a Gaussian probability distribution of these transformed 
variables appears to be more realistic than one for the correlation function itself. This expectation is verified with 
some numerical experiments which are described in Sect.lH Furthermore, we show that for a fiducial cosmological 
survey the Gaussian likelihood for the correlation function can significantly overlap with the region forbidden by ([2|), 
depending on survey size and number of separations at which the correlation function is measured. We conclude with 
a discussion and an outlook to future work and open questions. 



2. Upper and lower bounds on correlation functions 

Consider an n-dimensional homogeneous and isotropic random field g{x), with vanishing expectation value {g{x)) = 0, 
and with power spectrum -P(|A;|) and correlation function 

^P(|fe|)exp(-ifc.a:) , (3) 

which depends only on the absolute value of x, due to the assumed isotropy. This relation immediately shows that 

- m < a^) < m , (4) 

owing to P{k) > and |exp (— ife ■ x)\ < 1. However, the lower bound in Q is not an optimal one for more than one 
dimension. In two dimensions, the integral over the polar angle can be carried out, yielding 

r°° Ah k 

6-d(x)= -^P{k)Mkx), (5) 

where Jq is the Bessel function of the first kind of zero order. Since Jo(a;) has an absolute minimum at x sa 3.83 
with Jo, mill ~ —0.4028 (see also Abrahamsen 1997), the non-negativity of P(fc) impfies that ^2~d{x) > Jo,minC2-D(0). 
Similarly, in three dimensions one has 

e3-D(x)= / ^Pik)ioikx), (6) 
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where jo(a;) — sinx/x is the spherical Bessel function of zero order. Since jo (a;) has an absolute minimum at x w 4.493 
of jomin ~ -0.2172, the non-negativity of P{k) implies that ^3-0(2;) > jo,min 6-0(0). 
In the following, we will concentrate mainly on the one-dimensional case and write 

/•OO 

^{x) = dk Po{k) cos{xk) . (7) 
^0 

However, higher dimensions of the random field are included in all what follows, since by specifying x — (x,0, . . . ,0) 
in ([3]), we find 

P (fc) exp {-ikix) = J dfci cos(A:ia;) j^:^ J • ■ • dfc„ P{ki,k2,..., fc„) , (8) 

which thus takes the same form as ([7]). Thus, the n-dimensional case can be included in the same formalism as the 
one-dimensional case; note that for this argument, the random field is not restricted to be isotropic. However, as we 
shall discuss later, the resulting inequalities will not be optimal for isotropic fields of higher dimension. 

In the foregoing equations, P{k) can correspond either to the power spectrum of the underlying random process, or 
the sum of the underlying process and statistical noise. Furthermore, the power spectrum can also be the square of the 
Fourier transform of the realization of a random process in a finite sample volume. In all these cases, the non-negativity 
of P{k) applies, and the constraints on the corresponding correlation function derived below must hold. 



2.1. Constraints from the Cauchy-Schwarz inequality 

Making use of the Cauchy-Schwarz inequality. 



< / dfc f{k) / dfc h^{k) , (9) 



dfc f{k)h{k) 

we obtain by setting /(fc) = Po{k) and h{k) — ^ Po{k) co8{xk) thai0 

e{^) < m r dfc ^'o(fc) cos2(xfc) = m H ^fc ^'o(fc) ^^"""^^""^^ = ^ m + m^)] , m 



where we made use of the identity cos^ a = [1 + cos(2a)] /2. Together with ^ we therefore obtain the constraint 
equation 

-m + ^j^<m^)<m- (11) 

The interpretation of this constraint can be better seen in terms of the correlation coefficient r„ = ^(?ix)/^(0), which 
is defined for an arbitrary x. Then, (fTT|) reads 

-l + 2rl<r2<l. (12) 

This result can be interpreted as follows: If two points separated by x are strongly correlated, 1 — ri <C 1, then the 
value of the field at a position 2x must equally be correlated with that at x, which implies that also the correlation 
between the point 2x and the origin must be large. If the field at x is strongly anticorrelated with that at the origin, 
1 -I- ri <C 1, than the field at 2x must be similarly anticorrelated with that at a:, implying a strong correlation between 
the point 2x and the origin. The smaller |ri|, the weaker is the constraint (fT2)) . 
Making use of the identity 

[cos a 4- cos(2a)]"' = [1 + cos a] [1 + cos(3a)] , 

and applying the Cauchy-Schwarz inequality with f{k) = ^y Po{k)y/l + coa{xk) and h{k) — \/Poik)y^l + cos(3a;fc) , 
we find that 

[e(x) + a'2x)f < m) + a^)] im + a^^)] ■ (13) 

A second inequality is obtained by using in a similar way the identity 

[cos a — cos(2a)]^ = [1 — cos a] [1 — cos(3a)] . 
Both of these inequalities are summarized in terms of the correlation coefficient as 

(1 + ri) - (i-n) ^ ^ 



^ Note that this choice is possible because the power spectrum is non-negative! 
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Further inequahties involving ^(mx), with m > 4 being an integer, can be derived in this way. Making use of the 
relations 

[cos a + cos(n — l)a]^ [1 + cos(na)] [1 + cos(n — 2)a] ; [cos a — cos(n — l)a]^ — [1 ^ cos(na)] [1 — cos(n — 2)a] 
for n > 2, and employing the Cauchy-Schwarz inequality in the same way as before, we find 

- 1 + — ;— <rn<l , 15) 

1 + rn-2 1 - r„-2 

where the special case n = 3 has been derived already. We have thus found a set of inequalities for all correlation 
coefficients r„. In the next section, we will obtain bounds on the correlation function using a different method, and 
will show that these ones are stricter than those in (flSl) . 



2.2. Constraints from a covariance matrix approach 

We will proceed in a different way which is more straightforward. Consider a set of N points with m integer 

and < m < — 1. The covariance matrix of the random field at these N points has the simple form 

Cy = (g{ix)g{jx)) = - j\x) . (16) 

As is well known, the covariance matrix must be positive semi-definite, i.e., its eigenvalues must be non- negative. 
Dividing C by ^(0) > 0, we define 

A,, - Q,/m = n-jl - (17) 
and the eigenvalues of A must obey Xi > 0. For N = 2, the eigenvalues read Ai,2 = l±ri, yielding 

kil<l, (18) 

i.e. we reobtain For = 3, the eigenvalues read 



Ai — 1 — r2 , A2 3 — ^ 



2 

and the conditions Aj > can be solved for r2, yielding (H^). The four eigenvalues of A in the case iV = 4 are 

Al,2,3,4 ~ ^-^ 2 



ri + rs ± \ 5rf - 8rir2 + 4r2 - 2rir3 



and the conditions A; > after some algebra can be brought into the form For N > 5, the eigenvalues of A have 
a more complicated form; they are obtained as solutions of higher-order polynomials (see below) . 

However, we do not need an explicit expression for the eigenvalues, but only need to assure that they are non- 
negative. This condition can be formulated in a different way. The eigenvalues of the matrix A are given by the roots 
of the characteristic polynomial, which is the determinant of the matrix XSij — Aij. For a given N, this polynomial is 
of order in A and of the form 

JV-l 

A^ + y bkX" . (19) 



The coefficients bk of the polynomial are functions of the r^, as obtained from calculating the determinant. On the 
other hand, they can be expressed by the roots Xk of the polynomial; for example, for = 3, one finds that 

2 

A^ -f y bkX'' = X^ ~ (Ai + A2 + X3)X^ + (A1A2 + A1A3 + A2A3)A - A1A2A3 . 

The condition that all A/c > is then equivalent to the condition that 62 < 0, 61 > 0, and bo < 0. It is easy to show 
that these conditions lead to the constraint p^ . together with |ri| < 1. 

This procedure can be generalized for arbitrary A^: the condition that all eigenvalues of the correlation matrix A are 
non- negative is equivalent to requiring that the coefficients of the characteristic polynomial p9p satisfy b^-n < if n 
is odd, and b^-n > if is even. The explicit calculation of the characteristic polynomial by hand becomes infeasible, 
hence we use the computer algebra program Mathematica (Wolfram 1996). For successively larger A^, we calculate the 
characteristic polynomial. As we will show below, the characteristic polynomial factorizes into two polynomials; if A^ 
is even, it factorizes into two polynomials of order A^/2, whereas if A^ is odd, it factorizes into polynomials of order 
(A'' ± l)/2. The condition that all eigenvalues are non- negative is equivalent to the condition that the roots of both 
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polynomials are non-negative; this condition can be translated into inequalities for the coefficients of the polynomials 
in the same way as described above. 

We illustrate this procedure for the case N ^ 5. The characteristic polynomial in this case becomes 

(A^ + 6iA + 60) (A^ + caA^ + ciA + cq) , (20) 

with 

bi = r2 + r4 - 2 

bo = (l-r2)(l-r4)-(ri-r3)' 
C2 = ~3 - r2 - r4 

ci = 3(1 - rf) - 2rir3 - + 2r4 + r2(2 + n) - 2rl 

Co = {2rl - 1 - r2)r4 + rj + 2rir3(l - 2rl) + 2rl{l + r-i) - r2 + r?(3 - 4r2) - 1 

The conditions &i < 0, C2 < and ci > are satisfied, irrespective of the value of r4, owing to |r„| < 1 for all n, 
and the inequalities obtained before for r„, n < 3. Thus, these three inequalities do not yield additional constraints 
of r4. Those come from requiring 5o > and cq < 0. Both coefficients are linear in r4, which allows us to write the 
conditions explicitly, 

, , rf(l-4r2) + 2ri(l + r2) + 2rir3(l-2r2) + ri ^ ^ , (^^i - ^3)' 
+ l-2rl + r.2 (l-r2) ^''^ 

Alternatively, one can use the function InequalitySolve of Mathematica for the five inequalities for the four coeffi- 
cients, to obtain the same result. 

We see that the upper bound in (|2T|) agrees with that obtained from the Cauchy-Schwarz inequality - see (fTSj) 
- but the lower bounds differ. Hence, for n — A the bounds on r„ derived from the two methods are different, and 
that will be the case for all n > 4. This is not surprising, since the Cauchy-Schwarz inequality does not make any 
statement on how 'sharp' the bounds are, i.e., whether the bounds can actually be approached. One can argue, using 
discrete Fourier transforms, that the bounds obtained from the covariance method are 'sharp' (for a one-dimensional 
field) in the sense that for each combination of allowed r„, one can find a non-negative power spectrum corresponding 
to these correlation coefficients - the bounds can therefore not be narrowed further. 

It should be noted that in the cases given above, the lower (upper) bound on r„ is a quadratic function in r„_i, and 
its quadratic term cr'^_i has a positive (negative) coefficient c. This implies that the allowed region in the r„_i-r„-plane 
is convex; indeed, as we will show below, the allowed n-dimensional region for the ri, . . . , r„ is convex. 

The procedure illustrated for the case N — 5 holds for larger N as well: only the coefficients of A'' in the two 
polynomials yield new constraints on the correlation coefficient r^^i, and these two coefficients are linear in r^^i, 
so the constraints for r^-i can be obtained explicitly. This then leads us to conclude that the positive-definiteness of 
the matrix A for a given N is equivalent to the positivity of the determinants of all submatrices which are obtained 
by considering only the first n rows and columns of A, with 1 < n < A^. We will make use of this property in the next 
subsection. 

For larger A^, the upper and lower bounds are given by quite complicated expressions, so we refrain from writing 
them down; however, using the output from Mathematica, they can be used for subsequent numerical calculations. A 
much more convenient method for calculating the upper and lower bounds explicitly will be obtained below. 

To summarize, the procedure outlined gives us constraints on the form 

fnl < Tn < r„u , (22) 

where the lower and upper bounds are function of the r^, k < n, and satisfy —1 < r„i < r„u < 1. For n < 4, the 
bounds r„i and r„u have been written down explicitly above. 

The existence of these upper and lower bounds on r„ , and thus on the correlation function, immediately implies that 
the likelihood function of the correlation function cannot be a Gaussian, since the support of the latter is unbounded. 
This does not imply necessarily that the Gaussian approximation is bad; if the Gaussian probability ellipsoid described 
by 111) is 'small' in the sense that it is well contained inside the allowed region, the existence of the bounds alone yields 
no estimate for the accuracy of the Gaussian approximation. We will come back to this issue in Sect. [51 

Whereas the expressions for r„i and r„u for larger n become quite long, we found a remarkable result for the 
difference A„ = r„u ~ fni ■ This is most easily expressed by defining the functions 



n (^1 ; ■ ■ ■ ; ) — (^nu ^n) {j'n ) 



(23) 
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Then, for odd 

'(n-l)/2 \ /(«-l)/2 

A„(ri,...,r„_i) = 2 I -^2^ [] ^^k-i \ , (24) 



k=l / \ k=l 



and for even n, 

/ n/2 \ /n/2-1 

A„(ri,...,r„_i) = 2 ]^^2fe-i n -^2^1 (25) 



, fc=l / \ k=l 



This result has been obtained by guessing, and subsequently verified with Mathematica, up to n = 16, using the 
explicit expressions for the bounds derived next. We will make use of these properties in the Sect. [31 



2.3. Explicit expression for the bounds 

We will now derive an explicit expression for the upper and lower bounds on the r„. In doing so, we will first show 
that the determinant of the matrix A for N points factorizes into two polynomials, both of which are linear in r^-i- 
Then we make use of the fact that the positive definiteness of A is equivalent to having positive determinant of all 
submatrices obtained from A by considering only the first n rows and columns, n < N; however, these submatrices of 
A correspond to the matrix A for lower numbers of points, and their positive determinant is equivalent to the upper 
and lower bounds on the r„ for ti < — 2. Hence, by increasing N successively, the constraints on r^v-i are obtained 
from requiring det(A) > 0. 

The determinant of a matrix is unchanged if a multiple of one row (column) is added to another row (column). We 
make use of this fact for the calculation of det A, by carrying out the following four steps: 

1. We add the first row to the last, obtaining the matrix A^^' with elements 



2. We subtract the last column from the first. 



^(2) _ ^(1) _ . ^(1) 



The second row is added to the {N — l)-st one, the third row is added to the {N — 2)-nd one, and so on. For N 
odd, this reads 

(W-3)/2 



whereas for A'' even, the sum extends to {N — 2)/2. 
4. Finally, the {N — l)-st column is subtracted from the second one, the {N — 2)-nd column is subtracted from the 
third one, and so on, which for odd TV reads 

(7V-3)/2 

a(4) _ a(3) r a(3) 

k=l 

whereas for even N the sum extends to [N — 2)/2. 

These steps are illustrated for the case iV = 7 in Fig.[TJ For N even [odd], this results in a matrix A*^**^ which 
contains in the lower left corner a N/2 x N/2 [{N — l)/2 x {N + l)/2] submatrix with elements zero. Therefore, the 
determinant of A'^'*^ factorizes into the determinants of two N/2 x N/2 matrices [of a [N — l)/2 x [N — l)/2 and of a 
{N + l)/2 X (iV + l)/2 matrix]. Thus, 

det(A) = det (A^) = det (a^"*)) = det (U^) det (V^) , (26) 

where we explicitly indicate the number N of points, and thus the dimensionality of A, with a superscript, and where 
for even iV, and are N/2 x N/2 matrices with elements 

Ui} - r\^-J\ - rN+l-^-, ; V,^ = r|,_,| + r,+j_i , 1 < j < iV/2 , (27) 

whereas for odd N, and are (A^ - l)/2 x (iV - l)/2 and {N + l)/2 x {N + l)/2 matrices, respectively, with 
elements 

UiJ = r|,_,| - rAr+i_,_j , l<^,J<(^f-l)/2; V^} = r|,_,| + (1 - <5h) r,+,-2 , 1 < j < (A^ + l)/2 . (28) 
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A = 
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Fig. 1. For N = 7, the original covariance matrix A is shown, together with four transformations of it, described in 
the text, which leave the determinant unchanged 



Since Uf^ = l-rw-i, and the (7V/2-1) x (iV/2-1) (for TV even; for N odd this is a [(iV-l)/2-l] x [(TV- l)/2- 1]) 
submatrix obtained by cancelling the first column and row of is just U^^^, we can write 

det (U^) = (1 - rjv-i) det (U^-^) ^ (qn) ^ (29) 

where 0^ is the matrix which is obtained from by setting Uf^ = 0. The upper bound for rjv-i is found by setting 
det (U^) > 0, which yields 

det (U"+i) , , 

'^»<^"u = l + ^^^^, (30) 

where we set n = iV — 1. Analogously, the final element of V''^ reads = 1 + rjv-i, where m = N/2 for even N, 
and m={N + l)/2 for odd N. Therefore, 

det (V^) = (1 + rAT.i) det (y^-^) + det (V^) , (31) 

where is obtained from V-*^ by setting Vmm = 0; the lower bound for rjv-i is then obtained by setting this 
expression to zero, or 

det (V"+i) , ^ 

^"'=-^-d^tW- ^''^ 

Since det (U^) is a linear function of rjv-i, it must be of the form det (U^) = c(rjv-i — d), where the coefficients 
c,d are independent of tat-i- The value of d yields the root of det (U^), and is thus the upper bound on rN-i- The 
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coefficient c follows from therefore, 

dct (U"+i) = (r„u - r„) det (U"-i) . (33) 
The analogous result holds for the V^, i.e. 

det (V"+i) = (r„ - r„i) det (V"-i) . (34) 
These recursion relations then yield the explicit expressions 

7V/2 Af/2 

det (U^) - n (^(2fe-i)u - ^2fe-r) ; det (V^) - (^^fe-i - r(2fc-i)i) (35) 



fe=r k=i 



for TV even, and 



(7V-l)/2 (Af-l)/2 

det (U^) = n (r(2fc)u - r2fc) ; det (V^) = [] (ra^ - ^2k)i) (36) 

*;=! fe=l 

for N odd. This yields for the determinant of the matrix the explicit expression 

N/2 (JV-l)/2 

det(A^) =^.F2fe_i , det(A~)= ^ T^k (37) 

fe = l /£=1 

for even and odd N , respectively. Accordingly, the width A„ ~ r„u — r„i of the allowed range of r„ then becomes 

det(A^) , , 

where we made use of ([M)) and (^5]) . 
3. How strong are the constraints? 

We shall now investigate the question how constraining the constraints on the correlation function or, equivalently, the 
correlation coefficients are. In order to quantify this question, we imagine that the correlation function is measured on 
a regular grid of separations nx, and from that we define as before r„ = £,{nx)/^{0). Due to |r„| < 1 for all n. We 
therefore consider the set of functions defined by their values on the grid points, and allow only values in the interval 
— 1 < ''n < 1 for all n. Let M be the largest value of n that we consider; then the functions we consider are defined by 
a point in the M-dimensional space spanned by the r„ . This A/-dimensional cube has sidelength 2 and thus a volume 
of 2*^. We will now investigate which fraction of this volume corresponds to correlation functions which satisfy the 
constraints derived in the previous section. 

A related question that can be answered is: suppose the values of r^, 1 < A; < are given; what fraction of the 
(M — n)-dimensional subspace, spanned by the Vk with n + 1 < k < M corresponds to allowed correlation functions? 
For example, if ri = 1, then all other = 1, and hence the condition ri = 1 is very constraining - the volume fraction 
of the subspace spanned by the with A: > 2 is zero in this case. 

Mathematically, we define these volume fractions as 



2 / 2 ' ' ' / 2 

drfc 



n / ^- w 



k=n+l '''=1 



The factor 1/2 in each integration accounts for the side- length of the (M — n)-dimensional cube, so that the Vmu are 
indeed fractional volumes. Vmu depends on the with fc < n; in particular, Vmo = Vm is the fractional volume which 
is allowed if all constraints are taken into account. From the definition of the Vmu it is obvious that the following 
recursion holds: 

Vmu = / -^VMin+l) ■ (40) 

We can therefore calculate the Vun iteratively, starting with 

I drj\/ _ Ajv/ _ Tm-i^M'-z ■ ■ ■ 
yM{M-i) - / — - -:p :p , l^ij 
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where we have skipped the integration limits; here and in the following, an integral over r„ always extends from r„i 
to r„u- Here we made use of (l24l) or ([25)) . depending on M. The dots indicate that factors are added until !Fi or J^2 is 
reached. 

It should be noted that the only dependence of Vm(m-i) oh ^a/-i is through the factor Tm-i- Therefore, the next 
recursion step reads 



M(M-2) 



drM-1 



1 3^M-?,J-. 



V, 



M(Af-l) 



2 Tm-2^M-A ■ 
1 ^M-3^M-5 ■ 



2 J-M-2^M-i 

^M-2^M-i 



= 22b(2,2) 



AL_iB(2,2) 

2 



(42) 



(43) 

For the next step we notice that the only dependence of Vm(m-2) on rM-2 is through the function Tm-2^ which 
therefore lets us write 



where B(x, y) — V(x) V(y)lV{x + y) is the beta-function, and we used the relation 

dx {h - xy\x - a)" = (6 - a)i+2" B(l + n, 1 + n) 



M(M-3) 



drM- 



= 2B(2,2) 



Vm{M-2) 
!Fm-A^M-6 



M-2 



2^B(2,2)B(3,3) 



where we made again use of ()43p . Based on these results, we can now obtain a general expression, 

— n—2 ■ 



V, 



M(M- 



_ r,n(n-l) 



nB(fc,fc) 



.k=2 



J'M-n-lJ'M-n-3 



which can be proved by induction. In particular, the Vm = Vmo are given as 

M 



Y[Bik,k) 



.k=2 



(44) 



(45) 



(46) 



The values of Vm up to M — 20 are shown in Fig. [51 As can be seen from the upper panel, the admissible volume very 
quickly decreases as M increases. For comparison, in the lower panel we show V^/*^, i.e. the typical diameter of the 
allowed region. 



4. Generalizations 

The considerations of the previous sections were restricted to correlation functions as measured at equidistant points 
Xn = nx. In some cases, however, it may be advantageous to drop this constraint, e.g., to consider the correlation 
function at logarithmically spaced points. Furthermore, as mentioned before, tighter constraints on the correlation 
function are expected to hold in the multi-dimensional case. We shall consider these aspects in this section, presenting 
another method for deriving constraints which yields the optimal constraints for arbitrarily spaced points Xn in one- 
or higher-dimensional fields. But before, we briefly consider the covariance method for higher dimensions. 



4.1. Higher-dimensional fields: the covariance method 

As was noted above, the inequalities for the correlation functions have been obtained with a one-dimensional random 
field in mind. Whereas we have shown that all the bounds on correlation functions are also valid for higher-dimensional 
fields, they are not assumed to be 'optimal' - the reason is that the equivalent one-dimensional power spectrum defined 
in ([7|) was assumed to be totally arbitrary, except from being non-negative, whereas for isotropic fields in higher 
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M 

Fig. 2. Upper panel: volume fraction Vm (Eg. |46|) as function of the number M of separations for which the correlation 
function is measured. Lower panel: 1^^/*^ as function of M . One sees that the typical linear dimension of the allowed 
region decreases with M, meaning that the strong decrease of Vm with M is not just an effect of the dimensionality 
of the volume considered 



dimensions, it will obey further constraint relations due to the {n— l)-dimensional integration in ([8]). A first indication 
that the bounds can be improved has been seen with the lower bounds on the ratio ^(a::)/^(0), which turned out to 
be larger for two and three dimensions than for a one-dimensional field. The foregoing constraints on the correlation 
function are re-obtained with the covariance matrix approach in higher dimensions if the set of points are placed 
equidistant along one direction. New (and stronger) constraints are expected from this method if the distribution of 
points makes use of these higher dimensions. 

As a first example, we consider a two-dimensional (or higher-dimensional) field and place three points in an 
equilateral triangle of side-length x. The separation between any pair of points is then x, and the covariance matrix 
of these three points, normalized by ^(0), then reads 

(1 ri ri \ 
n 1 n . 
ri n 1 J 

The eigenvalues of this matrix are A1.2 = 1 — ri, A3 1 + 2ri, and requiring their non-negativity leads to 

- 1/2 < ri - a^)/m < 1 (47) 

for all X. The lower bound is somewhat smaller than the one obtained earlier for two-dimensional random fields, 
^5 —0.403, but significantly larger than that obtained for the 1-D case, ^(a;)/^(0) > —1. Next we consider a 
set of three points forming a triangle of which two sides have length x, and the third side has length rjx, with < 77 < 2. 
The corresponding covariance matrix reads 

(1 ri n \ 
n 1 / 

and its eigenvalues are Ai = 1 — r^, A2,3 = |^2 + ± ^ 8r^ + r^^ /2; note that we used the notation — ^{rix)/^{0). 
Non-negativity of the eigenvalues leads to the constraints 

max (-1/2, 2r2 - 1) < < 1 for < 7/ < 2 , (48) 

where we also used (|47t . 

Finally, we consider a square of side-length x, for which the covariance matrix reads 

/I n n 

ri 1 ''I ' 
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with the eigenvalues Ai_2 = 1 ^ A3, 4 = 1 + r ^ ± 1r\. Their non-negativity, combined with (|47p . yields 

max(-l/2,|ri|-l)<r^<l, (49) 

which is a weaker constraint than (|48| for r/ — \/2- Thus we see that the choice of the geometrical configuration 
of points affects the resulting constraints on the correlation function. It is by no means clear how to find a set of 
configurations such as to obtain the 'optimal' constraints in two (or higher) dimensions. In fact, methods other than 
using the covariance matrix need to be considered for obtaining optimal constraints (see below). 

Finally, we consider the simplest case in three dimensions, namely a set a of four points arranged in a regular 
tetrahedron of sidelength x. The resulting covariance matrix is 

/ 1 rx rx rx\ 

\ri ri ri 1 / 

with eigenvalues Ai^2,3 = 1 — ri, A4 = 1 + 3ri, yielding 

~ 1/3 < ri < 1 , (50) 

a constraint stronger than the ones obtained for one and two dimensions, but falling short of the one derived from the 
global minimum of the spherical Bessel function, ri —0.217. 



4.2. Optimal constraints in the general case 



We now describe a general method for deriving constraints on correlation functions £,{xn) of homogeneous and isotropic 
random fields, allowing for arbitrary values of separations Xn and arbitrary dimensions. However, it should be said 
right at the beginning that we were unable to obtain the corresponding bounds on the correlation functions explicitly. 
We can write the general relation ^ in the form 



ax) 



dk P{k) u{xk) 



(51) 



where P{k) is non-negative and m(0) = 1. For a 1-dimensional field, P{k) = P{k)/ {2tt), u{y) = cosy; for two dimensions, 
P{k) = k P{k)/ {2it), u{y) — Jo(y), and for three dimensions, P{k) = fc^ P(fc)/(27r^), u{y) = jo(y). Next we consider a 
quadrature formula for the integral, and write 



K 



K 



(52) 



where the Wj are (positive) weights corresponding to the quadrature formula, and in the last step we defined Wj > O. 
This approximation can be made arbitrarily accurate by letting K ^ 00. Defining the correlation coefhcient as before, 
we obtain 

K 

r{x)^i{x)/m^Y.'^3<xk,), (53) 



where the coefficients 



K 



(54) 



satisfy O < < 1 and = 1- 

If we now consider a set of N points a;„, together with the correlation coefficients r„ = r(x„), then we see that 
a point in the iV-dimensional space of r = (ri, . . . , rjv) can be described as a weighted sum of points lying along the 
curve c(A) = (u(Axi), . . . , u{\xn)), < A < 00, i.e.. 



K 



(55) 



where we considered the transition of the discrete points kj to a continuous variable A. Since Vi G [0, 1] and ^ 1^ = 1, 
the point r must be located inside the convex volume containing all points on the curve c(A) . It is clear that this convex 
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Fig. 3. Left panel: The curve c(A) — (Jo(A), Jo(2A)), together with its convex envelope. The latter consists of two 
segments of the curve c and two straight lines, one of which is tangent to c(A) at the points c(Ai) sa (0.1905, —0.3870) 
and c(A2) ~ (-0.1894, -0.2421), the other one is tangent to c at c(A3) w (-0.3872, 0.2987) and intersects c(0) = (1, 1). 
The corresponding values of the curve parameter are Ai w 2.0580, A2 ~ 4.9636, and A3 ~ 3.5561. These numerical 
values are found as follows: the requirement that the lower straight line segment is tangent to curve at the two points 
Ai^2 leads to the conditions c(Ai) = [c(Ai) — c(A2)] for i = 1,2, where are scalars. These four scalar equations 
for the four unknowns a,;, A^ have several solutions, the relevant one is the 'outermost' one. For the upper straight 
line segment, one employs the condition that the tangent at point A3 goes through c(0), i.e., c(A3) = a [c(A3) — c(0)], 
and of the multiple solutions of these two scalar equations for the two unknowns a and A3, one takes the 'outermost' 
one. Middle panel: in a similar way, the curve c(A) = (jo(-^)) jo(2A)) is plotted, together with its convex envelope. 
Here, Ai « 2.3911, A2 ~ 5.3490 and A3 w 4.0287, and c(Ai) w (0.2852,-0.2086), c(A2) w (-0.1503,-0.0894), and 
c(A3) w (—0.1924, 0.1216). Right panel: Comparison of the allowed regions in the ri-r2-plane in 1, 2 and 3 dimensions. 



envelope of the curve c yields indeed the optimal general bounds on the r^: every point within the convex envelope 
can be realized by choosing a set of n points on the curve c appropriately. Since the function u(y) depends on the 
dimension of the random field, the constraints will be different for different numbers of dimensions. Furthermore, the 
curve c(A) depends on the choice of the x„, and therefore the constraints will also depend on the choice of separations 
for which the correlation function is measured. 

Unfortunately, we have not found a way how to algebraically describe the convex envelope of the curve c(A), and 
thus to obtain explicit expressions for the upper and lower bounds on r.i. For now, we therefore will present just a few 
simple examples. 

For the 1-dimensional case with two points = 2x1, the curve reads c(A) = (cos A, cos 2A), with A ~ x\k\ the 
same set of points is described by the curve c'(a) — {a,2a^ — 1) (using trigonometric identities), — 1 < a < 1. Thus, 
the convex envelope of c' is the region between the parabola c' and the line r2 — +1, and thus we re-obtain the 
bounds pT|) . For the choice X2 = 3a;i, the set of points of c(A) is equivalent to that of the curve c'(a) = (a, 4a^ — 3a), 
a € [—1, 1]. The convex envelope can then be described by the bounds max(— l,4rj — 3ri) < r2 < min(l,4rij' — 3ri), 
where the lower bound differs from —1 for ri > 1/2, and the upper bound is different from +1 for ri < —1/2. If 
we choose instead X2 — ^ixi, then for a generic (non-rational) /i, the curve c(A) fills the whole square —1 < ri < 1, 
— 1 < ''2 < Ij which is also coincident with its convex envelope. 

As the next example, we consider a 2-dimensional field with X2 — 2xi. In the left panel of Fig. [31 the curve 
c(A) = (Jo(A), Jo(2A)) is plotted for A > 0. The convex envelope in this case can be constructed explicitly: The 
boundary of the smallest convex region which contains the curve c is composed of four parts: (1) The section of the 
curve c(A) for < A < Ai, (2) the straight line connecting the two points c(Ai) and c(A2), (3) the section of the 
curve c(A) for A3 < A < A2, and (4) the straight line connecting c(A3) and c(0) — (1, 1). In a similar way, the convex 
envelope can be constructed for a 3-dimensional random field; see Fig. [31 middle panel. 

Unfortunately, we have not yet found a systematic way how to construct the convex envelope for n > 2, and how 
to obtain explicit bounds on the correlation coefficients in these cases - although it is clear that the allowed region, 
at least in the r'i-r2-plane, decreases as one goes to higher-dimensional fields (see right panel of Fig.[3|). Therefore, the 
development of explicit bounds in higher dimensions is of great interest. 
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5. Transformation of variables 

The finite bounds on the correlation coefficients clearly show that the likelihood of the correlation function cannot 
be a Gaussian. However, the bounds on r„ may suggest that the Gaussian approximation for the likelihood could be 
better in terms of transformed variables, in which the allowed range for each ?■„ is mapped onto the real axis. Such a 
transformation would simplify the parametrization of the likelihood from numerical simulations. Defining 

^71U '^n\ /r-r'\ 

Xn := , (56) 

the allowed range of r„ is mapped onto — 1 < < 1. It should be noted that (j56p is a coupled non- linear transformation 
of variables, since the bounds on r„ depend on the rk, k < n. The mapping to the real axis is then obtained by any 
function that maps [—1, 1] to (— cx),oo); we choose 

yn ■■= atanh(a;„) . (57) 

Thus, no bounds on the correlation coefficients are violated if the probability distribution of the ?/„ would follow a 
multi-variate Gaussian. 

We next consider the relation between the likelihood of the correlation function and the corresponding likelihood 
of the Un- To shorten notation, we drop the explicit dependence on the model parameters p in the following. Then, 

N N N 



fe=0 k=l k=l 

N 

= i^v{yi,---,yN\£.o)W_Ayk Co{£,o)d£,o , (58) 



fc=i 



where in the first step we have used the product rule of probability theory and introduced the conditional probability 
density C , and where the next two steps define the likelihood in terms of the Vn and the j/n, respectively. Thus, the 
distribution of the r„ and yn can depend explicitly on the value of ■ 



The relation between Cr and Cy is obtained from 



Lr{ri, . . . ,rAr|^o) = Cy{yi, . . . , yAr|^o) det( J) , (59) 
where the yi are functions of the rj , and the transformation matrix J is given by 



_ dyi 1 dx, 2 



~ dr, 1 - xj dr, Af - (2r, - - vaf 



(60) 



Note that the partial derivatives vanish for j > i. Thus, the transformation between the probability distribution of the 
Ti and that of the yi is rather complicated, implying that the behaviour of these two probability distributions can be 
considerably different. In particular, the Cy could be approximated by a Gaussian, the corresponding Cr would have 
a significantly different functional form; also the covariances of the two distributions will differ significantly. 



6. Simulations 

In order to illustrate the analytical results discussed in the previous sections and to explore the effects of the constraints 
on the shape of the likelihood of the correlation function, we have conducted some simple numerical experiments. We 
have generated realizations of a periodic one-dimensional Gaussian random field 5 on a regular grid according to 

1 ^ 

5'-^E«''^"'''^''5,' (61) 

where gi is the discrete Fourier transform of the random field, and has a normal distribution of zero mean, gi ~ 
Af (0,y/N where Pi is the discretized power spectrum and M{iJi,a) is the Gaussian distribution with mean /i 

and standard deviation a. From each realization, we measure the correlation function using the estimator 

1 ^ 

ii = J^"^ 9a ga+t (62) 

and calculate the correlation coefficients r^. Note that this estimator of the correlation function explicitly employs 
the periodicity of the discrete random field. As defined in this way, the estimated correlation function is the Fourier 
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Fig. 4. Constraints in the ri-r2-plane (left panel) and in the r2-r3-plane (left panel), where ri was constrained to 
—0.41 < ri < —0.39. Each point corresponds to a correlation function measured from a realization of a one-dimensional 
Gaussian random field with a randomly drawn power spectrum and = 16 modes. Plotted as lines are the analytically 
determined constraints as given by Eqs. (fT2ll and (fT4|) . 




Fig. 5. Distribution of (ri, r2) for correlation functions measured from simulated Gaussian random fields with N — 256 
and random power spectra (similar to Fig. [4]). Even though the points lie well inside the admissible range, the shape 
of the distribution is similar to the shape of the allowed region. 

transform of a power spectrum P/, proportional to the squares of the amplitudes of the gi. Hence, P- is non- negative, 
and we thus expect that the correlation function of each realization obeys the constraints derived before - indeed, this 
is the case. 

We illustrate the constraints in the ri-r2-plane (Eq. [T^ and in the r2-r3-plane for fixed ri (Eq. [T?|) in Fig. U) In 
order to fully populate the allowed regions with data points, we do not use a single power spectrum for all realizations 
of the random field, but draw for each realization random positive values, uniformly distributed in [0, 1], for each Pi. 
Note that the data points do not completely fill out the regions allowed by Eqs. (fT2|) and (fT4|) . but are constrained 
to a slightly smaller region with piecewise linear boundaries. The reason for this is that we have used a random field 
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with N = 16 modes, whereas Eqs. ([T^ and have been derived without reference to the specific way the random 
field is created and are vahd also for an infinite number of modes. 

The covariance matrix of the correlation function estimator given by Eq. (|62p is given by 

^ N N 

Cov[^]y = J^'^Yl (^\"'-''\ ^\a+t-b-3\ + C\a-b-3\ ^\a-b+t\) (63) 
a=l 6=1 

With this, we can compare the true distribution of the correlation function (as estimated from a large number of 
realizations of the Gaussian field) to the commonly used Gaussian approximation to the likelihood. As an example, we 
show the likelihoods /^(Csj^io) Siiid /^(CiOjCis) Fig. [SI where we have marginalized over all other components of ^. 
We have used a random field with A'^ = 64 modes and a single Gaussian power spectrum. Even though the estimated 
likelihood and the Gaussian prediction have identical first and second moments, a Gaussian likelihood clearly is a bad 
approximation to the distribution of the correlation function. 

We now investigate whether the transformation of variables from r„ to ?/„ given by Eqs. (|56p and (|57p leads to a 
likelihood of the y„ that is better described by a Gaussian. We find that this is indeed the case, as is illustrated in 
Figs. [7] and [5] for the marginalized distributions in the 1-2-, and 2-3-planes. In the left panel, we show the estimated 
distribution of the correlation coefhcients, whereas the right panel contains the distribution after the transformation. 
The probability distribution in r-space 'feels' the presence of the boundary between allowed and forbidden regions, 
even at the innermost contours which are quite well away from the boundary. We also show the contours of the best- 
fitting bi-variate Gaussian for comparison and see that the distribution in the y-space is much better approximated by 
a Gaussian than the distribution of the original correlation coefficients. In addition, we note that the transformation 
r t—f y seems to reduce the correlations between most of the components of the correlation function. It may be 
suggested that the fact of the far more Gaussian distribution in the transformed variable is related to the choice of 
the transformation whose functional form is that of Fisher's z-transformation (Fisher 1915). 

Finally, we wish to assess the importance of the constraints in a more realistic, two-dimensional setting and consider 
as an example a weak lensing survey. Our strategy is as follows: we draw a large number of realizations of the shear 
correlation function ^_|_, which is the two-dimensional Fourier transform of the convergence power spectrum (Kaiser 
1992; Bartelmann & Schneider 2001) from a multivariate Gaussian likelihood. The covariance matrix for the likelihood 
function is computed under the assumption that the convergence is a Gaussian random field using the methodology 
described in Joachimi et al. (2008). For each of these realizations, we compute the matrix A (see Eq. [T7| . We then use 
this sample of correlation functions to compute a Monte-Carlo estimate of the integral 

A = y d"C£(Cb)i7pd(A) , (64) 

where C{£,\p) denotes the Gaussian likelihood as given by Eq. ^ and ffpd(A) = 1 if A is positive semi-definite and 
-ffpd(A) = otherwise. We test A for positive semi-dcfinitcness using the Cholesky-decomposition (Press et al. 1992). 
A measures the overlap of the Gaussian distribution with the allowed region. If A < 1, the Gaussian likelihood assigns 
a finite probability to regions containing correlation functions that do not correspond to a positive power spectrum 
and are therefore forbidden by the constraints discussed in this paper. We can therefore use A as a rough indicator 
of the validity of the assumption of a Gaussian likelihood. However, note that even if A « 1, the shape of the true 
likelihood might deviate significantly from a Gaussian distribution. 

For the numerical experiment, we choose a WMAP-5-like cosmology to compute the shear correlation function and 
its covariance matrix (approximating the shear field to be Gaussian) , and a redshift distribution of the source galaxies 
similar to the one found for the CFHTLS-Wide (Benjamin et al. 2007). The results of the numerical experiment 
are shown in Fig. [9l where we plot A as a function of the number of bins n (linear binning), keeping the maximum 
angular scale to which ^_|- is evaluated constant at 6'inax = 20'. We display the results for three different survey sizes 
(1, 10 and 100 deg^). The constraints are more important for smaller survey areas and a larger number of bins. This 
is because as the number of bins increases, the constaints on the admissible values for the following components 
of the correlation function become tighter, so that the correlation function is basically determined by its first few 
components. Increasing the noise (small area) or increasing the number of bins therefore decreases the fraction of 
positive semi-definite correlation functions that can be drawn from the Gaussian likelihood. These results indicate 
that the likelihood of the correlation function might deviate significantly from a Gaussian even in realistic situations. 
It is expected that for real shear fields, which are not Gaussian on the angular scales considered here, the values of A 
will deviate from unity even more. We speculate that the non-Gaussianity found in Hartlap et al. (2009) for the shear 
correlation functions might at least partly be caused by the constraint of positive semi-definiteness. 
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Fig. 6. Marginalized likelihood of the correlation function components ^5 and ^10 (left panel) and ^10 and ^15 (right 
panel) as estimated from 10^ realizations of a one-dimensional Gaussian random field with iV = 64 modes and a 
Gaussian power spectrum (black contours). The corresponding Gaussian likelihood with covariance matrix as predicted 
from the power spectrum (Eq. [63]) is shown by the red dashed contours. 




ri yi 

Fig. 7. Marginalized distribution of ri and r2 for a Gaussian random field with A'^ = 64 modes and Gaussian power 
spectrum (right panel) and the distribution of the corresponding transformed variables yi and y2 (right panel) . Shown 
as red dashed contours are the best-fitting bi-variate Gaussian distributions; in the left panel, the constraint given by 
Eq. ([T2I) is shown as solid blue line. 



7. Conclusions and outlook 

We have considered constraints on correlation functions that need to be satisfied in order for the correlation function to 
correspond to a non-negative power spectrum. Using the covariance matrix method, we have derived explicit expressions 
for the upper and lower bounds on the correlation coefficients r„; these were derived for the case that the spatial 
sampling of occurs at points Xj = jx. This method yields optimal constraints for the correlation of one-dimensional 
random fields, whereas they are not optimal for higher-dimensional homogeneous and isotropic random fields. We 
have indicated a method with which such optimal constraints can in principle be derived for higher-dimensional fields 
and for non-linear spacing of grid points, and presented a few simple applications of this method; however, up to now 
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1*2 yi 

Fig. 8. Distribution of r2 and for a Gaussian random field with TV = 64 modes and Gaussian power spectrum (right 
panel), where we set 0.05 < ri < 0.15 and have marginalized over all other components of the correlation function, 
and the distribution of the corresponding transformed variables ?/2 and j/3 (right panel) . Shown as red dashed contours 
are the best-fitting bi-variate Gaussian distributions; in the left panel, the constraint given by Eq. for ri =0.1 is 
shown as solid blue line. 
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we have not obtained a systematic method how to derive explicit upper and lower bounds on the r„ in these cases. 
Finding those will be of considerable interest since they are expected to be tighter than the corresponding bounds 
derived for the 1-D case. 

Using cosmic shear as an example, we have demonstrated that the Gaussian probability ellipsoid, which is obtained 
under the assumption that the likelihood function of the correlation function is given by a Gaussian characterized by 
the covariance matrix, significantly spills over to the forbidden region of correlation functions. This effect is even more 
serious than considered here, since we have used for this analysis the constraints from Sect. [5] which are not optimal. 
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as shown in Sect.[31 Hence, from this argument alone we conclude that the assumption of a Gaussian likelihood is not 
very realistic and probably lead to erroneous estimates of parameters and their confidence regions. 

Even if the Gaussian likelihood ellipsoid is contained inside the allowed region, the true likelihood deviates from 
a Gaussian, as simple numerical experiments with one-dimensional Gaussian random fields have shown. Surprisingly, 
the shape of the resulting likelihood contours have some resemblance to the shape of the boundary of the allowed 
region even if the size of the probability distribution is considerably smaller than the allowed region. The origin of 
this result is not understood. Most likely, the likelihood of the correlation function depends not only on its covariance 
matrix, but also on higher-order moments. 

A non-linear coupled transformation of the correlation coefficients leads to a distribution that appears much more 
Gaussian (in the transformed variables), and there may be a connection of this fact to the Independent Components 
analyzed in Hartlap et al. (2009). Indeed, in these transformed variables, the likelihood not only is much more Gaussian, 
but also less correlated, which supports the hypothesis about a connection between the constraints derived here and 
the ICA analysis of Hartlap et al. (2009). More extensive numerical tests may yield better insight into this connection. 
It must be stressed that such a result, if it can be obtained, would be of great importance, given that the determination 
of multi-variate probability distributions from numerical simulations is prohibitively expensive. 

An alternative route for understanding the connection between the constraints derived here and the shape of the 
likelihood function is the explicit calculation of the multivariate probability distribution of the correlation function 
£,{xj) for a Gaussian random field. The constraints on the rj should be explicitly present in this probability distribution. 
Work on these issues is currently in progress. 

The results of this paper can most likely be generalized to random fields which are not scalars, e.g., the polarization 
of radiation, or the orientation of objects. The cosmic shear correlation function (e.g.. Kaiser 1992; Schneider et al. 
2002) ^_|_ which has been considered in Sect. [S] is equivalent to the correlation function of the underlying surface mass 
density k, and thus the correlation behaves in the same way as that of a scalar field. However, the other cosmic shear 
correlation function ^_ is qualitatively different, being a spin-4 quantity, and for which the filter function in ^ is 
replaced by Ji{kx). 

The foremost aim of this paper was the derivation of exact constraints on correlation functions; in contrast, we have 
not considered methods for measuring a correlation function from data. For example, in many cases the correlation 
function cannot be measured at zero lag, so that the correlation coeffcients r = £,{x)/^(0) can then not be determined 
directly. Furthermore, one derives the correlation function from data in a given volume, and thus the measured 
correlation function will deviate from the ensemble average, even in the absence of noise. This effect has two different 
aspects: suppose for a moment that the observed field is one-dimensional and forced (or assumed) to be periodic. If the 
correlation function on such a field is measured using the definition (p^ . then the measured correlation function will 
deviate from the ensemble average, but it will still correspond to a non-negative power spectrum, given by the square 
of the Fourier transform of the realization of the field. Hence, in such a case, every measured correlation function will 
satisfy the constraints derived in Sect. [21 If the field has more than one dimension, but is still periodic, the correlation 
function measured by a generalization of ([62]) to higher dimensions will still obey the constraints from Sect. [21 for the 
same reason. However, the power spectrum of the realization of the field will in general not be isotropic, and thus the 
considerations of Sect.Hl do not apply strictly. In the more realistic case where periodicity cannot be assumed, one 
cannot measure the correlation function by 'wrapping around' as in l|62p : in this case, there are 'boundary effects', 
which are the stronger the more the separation approaches the size of the data field. Then, there exists not necessarily 
a non-negative power spectrum related to the measured correlation function through and the constraints may not 
apply strictly. How important these effect are needs to be studied with a more extended set of simulations. 
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